home *** CD-ROM | disk | FTP | other *** search
/ Programmer Power Tools / Programmer Power Tools.iso / math / praxis.arc / PRAXIS.C < prev    next >
C/C++ Source or Header  |  1987-07-15  |  17KB  |  540 lines
  1. /*********************************************************************/
  2. /*     f u n c t i o n     p r a x i s                              */
  3. /*                                                                   */
  4. /* praxis is a general purpose routine for the minimization of a     */
  5. /* function in several variables. the algorithm used is a modifi-    */
  6. /* cation of conjugate gradient search method by powell. the changes */
  7. /* are due to r.p. brent, who gives an algol-w program, which served */
  8. /* as a basis for this function.                                     */
  9. /*                                                                   */
  10. /* references:                                                       */
  11. /*     - powell, m.j.d., 1964. an efficient method for finding       */
  12. /*       the minimum of a function in several variables without      */
  13. /*       calculating derivatives, computer journal, 7, 155-162       */
  14. /*     - brent, r.p., 1973. algorithms for minimization without      */
  15. /*       derivatives, prentice hall, englewood cliffs.               */
  16. /*                                                                   */
  17. /*     problems, suggestions or improvements are always wellcome     */
  18. /*                       karl gegenfurtner   07/08/87                */
  19. /*                                           c - version             */
  20. /*********************************************************************/
  21. /*                                                                   */
  22. /* usage: min = praxis(fun, x, n);                                   */
  23. /*                                                                   */
  24. /*  fun        the function to be minimized. fun is called from      */
  25. /*             praxis with x and n as arguments                      */
  26. /*  x          a double array containing the initial guesses for     */
  27. /*             the minimum, which will contain the solution on       */
  28. /*             return                                                */
  29. /*  n          an integer specifying the number of unknown           */
  30. /*             parameters                                            */
  31. /*  min        praxis returns the least calculated value of fun      */
  32. /*                                                                   */
  33. /* some additional global variables control some more aspects of     */
  34. /* the inner workings of praxis. setting them is optional, they      */
  35. /* are all set to some reasonable default values given below.        */
  36. /*                                                                   */
  37. /*   prin      controls the printed output from the routine.         */
  38. /*             0 -> no output                                        */
  39. /*             1 -> print only starting and final values             */
  40. /*             2 -> detailed map of the minimization process         */
  41. /*             3 -> print also eigenvalues and vectors of the        */
  42. /*                  search directions                                */
  43. /*             the default value is 1                                */
  44. /*  tol        is the tolerance allowed for the precision of the     */
  45. /*             solution. praxis returns if the criterion             */
  46. /*             2 * ||x[k]-x[k-1]|| <= sqrt(macheps) * ||x[k]|| + tol */
  47. /*             is fulfilled more than ktm times.                     */
  48. /*             the default value depends on the machine precision    */
  49. /*  ktm        see just above. default is 1, and a value of 4 leads  */
  50. /*             to a very(!) cautious stopping criterion.             */
  51. /*  step       is a steplength parameter and should be set equal     */
  52. /*             to the expected distance from the solution.           */
  53. /*             exceptionally small or large values of step lead to   */
  54. /*             slower convergence on the first few iterations        */
  55. /*             the default value for step is 1.0                     */
  56. /*  scbd       is a scaling parameter. 1.0 is the default and        */
  57. /*             indicates no scaling. if the scales for the different */
  58. /*             parameters are very different, scbd should be set to  */
  59. /*             a value of about 10.0.                                */
  60. /*  illc       should be set to true (1) if the problem is known to  */
  61. /*             be ill-conditioned. the default is false (0). this    */
  62. /*             variable is automatically set, when praxis finds      */
  63. /*             the problem to be ill-conditioned during iterations.  */
  64. /*  maxfun     is the maximum number of calls to fun allowed. praxis */
  65. /*             will return after maxfun calls to fun even when the   */
  66. /*             minimum is not yet found. the default value of 0      */
  67. /*             indicates no limit on the number of calls.            */
  68. /*             this return condition is only checked every n         */
  69. /*             iterations.                                           */
  70. /*                                                                   */
  71. /*********************************************************************/
  72.  
  73. #include <math.h>
  74. #include <stdio.h>
  75. #include "machine.h"
  76.  
  77.  
  78. /* control parameters */
  79. double tol = SQREPSILON,
  80.        scbd = 1.0,
  81.        step = 1.0;
  82. int    ktm = 1,
  83.        prin = 2,
  84.        maxfun = 0,
  85.        illc = 0;
  86.        
  87. /* some global variables */
  88. static int i, j, k, k2, nl, nf, kl, kt;
  89. static double s, sl, dn, dmin,
  90.        fx, f1, lds, ldt, sf, df,
  91.        qf1, qd0, qd1, qa, qb, qc,
  92.        m2, m4, small, vsmall, large, 
  93.        vlarge, ldfac, t2;
  94. static double d[N], y[N], z[N],
  95.        q0[N], q1[N], v[N][N];
  96.  
  97. /* these will be set by praxis to point to it's arguments */
  98. static int n;
  99. double *x;
  100. double (*fun)();
  101.  
  102. /* these will be set by praxis to the global control parameters */
  103. static double h, macheps, t;
  104.  
  105. double 
  106. random()    /* return random no between 0 and 1 */
  107. {
  108.    return (double)(rand()%(8192*2))/(double)(8192*2);
  109. }
  110.  
  111. sort()        /* d and v in descending order */
  112. {
  113.    int k, i, j;
  114.    double s;
  115.  
  116.    for (i=0; i<n-1; i++) {
  117.        k = i; s = d[i];
  118.        for (j=i+1; j<n; j++) {
  119.            if (d[j] > s) {
  120.           k = j;
  121.           s = d[j];
  122.        }
  123.        }
  124.        if (k > i) {
  125.       d[k] = d[i];
  126.       d[i] = s;
  127.       for (j=0; j<n; j++) {
  128.           s = v[j][i];
  129.           v[j][i] = v[j][k];
  130.           v[j][k] = s;
  131.       }
  132.        }
  133.    }
  134. }
  135.  
  136. print()        /* print a line of traces */
  137. {
  138.    int i;
  139.  
  140.    printf("\n");
  141.    printf("... chi square reduced to ... %20.10e\n", fx);
  142.    printf("... after %u function calls ...\n", nf);
  143.    printf("... including %u linear searches ...\n", nl);
  144.    vecprint("... current values of x ...", x, n);
  145. }
  146.  
  147. matprint(s, v, n)
  148. char *s;
  149. double v[N][N];
  150. {
  151.    int k, i;
  152.    
  153.    printf("%s\n", s);
  154.    for (k=0; k<n; k++) {
  155.        for (i=0; i<n; i++) {
  156.            printf("%20.10e ", v[k][i]);
  157.        }
  158.        printf("\n");
  159.    }
  160. }
  161.  
  162. vecprint(s, x, n)
  163. char *s;
  164. double x[N];
  165. {
  166.    int i;
  167.    
  168.    printf("%s\n", s);
  169.    for (i=0; i<n; i++)
  170.        printf("%20.10e ", x[i]);
  171.    printf("\n");
  172. }
  173.  
  174. #ifdef MSDOS
  175. static double tflin[N];
  176. #endif
  177.  
  178. double
  179. flin(l, j)
  180. double l;
  181. {
  182.    int i;
  183. #ifndef MSDOS
  184.    double tflin[N];
  185. #endif   
  186.  
  187.    if (j != -1) {        /* linear search */
  188.       for (i=0; i<n; i++)
  189.           tflin[i] = x[i] + l *v[i][j];
  190.    }
  191.    else {            /* search along parabolic space curve */
  192.       qa = l*(l-qd1)/(qd0*(qd0+qd1));
  193.       qb = (l+qd0)*(qd1-l)/(qd0*qd1);
  194.       qc = l*(l+qd0)/(qd1*(qd0+qd1));
  195.       for (i=0; i<n; i++)
  196.           tflin[i] = qa*q0[i]+qb*x[i]+qc*q1[i];
  197.    }
  198.    nf++;
  199.    return (*fun)(tflin, n);
  200. }
  201.  
  202. min(j, nits, d2, x1, f1, fk)
  203. double *d2, *x1, f1;
  204. {
  205.    int k, i, dz;
  206.    double x2, xm, f0, f2, fm, d1, t2,
  207.           s, sf1, sx1;
  208.  
  209.    sf1 = f1; sx1 = *x1;
  210.    k = 0; xm = 0.0; fm = f0 = fx; dz = *d2 < macheps;
  211.    /* find step size */
  212.    s = 0;
  213.    for (i=0; i<n; i++) s += x[i]*x[i];
  214.    s = sqrt(s);
  215.    if (dz)
  216.       t2 = m4*sqrt(fabs(fx)/dmin + s*ldt) + m2*ldt;
  217.    else
  218.       t2 = m4*sqrt(fabs(fx)/(*d2) + s*ldt) + m2*ldt;
  219.    s = s*m4 + t;
  220.    if (dz && t2 > s) t2 = s;
  221.    if (t2 < small) t2 = small;
  222.    if (t2 > 0.01*h) t2 = 0.01 * h;
  223.    if (fk && f1 <= fm) {
  224.       xm = *x1;
  225.       fm = f1;
  226.    }
  227.    if (!fk || fabs(*x1) < t2) {
  228.       *x1 = (*x1 > 0 ? t2 : -t2);
  229.       f1 = flin(*x1, j);
  230.    }
  231.    if (f1 <= fm) {
  232.       xm = *x1;
  233.       fm = f1;
  234.    }
  235. next:
  236.    if (dz) {
  237.       x2 = (f0 < f1 ? -(*x1) : 2*(*x1));
  238.       f2 = flin(x2, j);
  239.       if (f2 <= fm) {
  240.          xm = x2;
  241.      fm = f2;
  242.       }
  243.       *d2 = (x2*(f1-f0) - (*x1)*(f2-f0))/((*x1)*x2*((*x1)-x2));
  244.    }
  245.    d1 = (f1-f0)/(*x1) - *x1**d2; dz = 1;
  246.    if (*d2 <= small) {
  247.       x2 = (d1 < 0 ? h : -h);
  248.    }
  249.    else {
  250.       x2 = - 0.5*d1/(*d2);
  251.    }
  252.    if (fabs(x2) > h)
  253.       x2 = (x2 > 0 ? h : -h);
  254. try:
  255.    f2 = flin(x2, j);
  256.    if ((k < nits) && (f2 > f0)) {
  257.       k++;
  258.       if ((f0 < f1) && (*x1*x2 > 0.0))
  259.          goto next;
  260.       x2 *= 0.5;
  261.       goto try;
  262.    }
  263.    nl++;
  264.    if (f2 > fm) x2 = xm; else fm = f2;
  265.    if (fabs(x2*(x2-*x1)) > small) {
  266.       *d2 = (x2*(f1-f0) - *x1*(fm-f0))/(*x1*x2*(*x1-x2));
  267.    }
  268.    else {
  269.       if (k > 0) *d2 = 0;
  270.    }
  271.    if (*d2 <= small) *d2 = small;
  272.    *x1 = x2; fx = fm;
  273.    if (sf1 < fx) {
  274.       fx = sf1;
  275.       *x1 = sx1;
  276.    }
  277.    if (j != -1)
  278.       for (i=0; i<n; i++)
  279.           x[i] += (*x1)*v[i][j];
  280. }
  281.  
  282. quad()    /* look for a minimum along the curve q0, q1, q2    */
  283. {
  284.    int i;
  285.    double l, s;
  286.  
  287.    s = fx; fx = qf1; qf1 = s; qd1 = 0.0;
  288.    for (i=0; i<n; i++) {
  289.        s = x[i]; l = q1[i]; x[i] = l; q1[i] = s;
  290.        qd1 = qd1 + (s-l)*(s-l);
  291.    }
  292.    s = 0.0; qd1 = sqrt(qd1); l = qd1;
  293.    if (qd0>0.0 && qd1>0.0 &&nl>=3*n*n) {
  294.       min(-1, 2, &s, &l, qf1, 1);
  295.       qa = l*(l-qd1)/(qd0*(qd0+qd1));
  296.       qb = (l+qd0)*(qd1-l)/(qd0*qd1);
  297.       qc = l*(l+qd0)/(qd1*(qd0+qd1));
  298.    }
  299.    else {
  300.       fx = qf1; qa = qb = 0.0; qc = 1.0;
  301.    }
  302.    qd0 = qd1;
  303.    for (i=0; i<n; i++) {
  304.        s = q0[i]; q0[i] = x[i];
  305.        x[i] = qa*s + qb*x[i] + qc*q1[i];
  306.    }
  307. }
  308.  
  309. double
  310. praxis(_fun, _x, _n)
  311. double (*_fun)();
  312. double _x[N];
  313. {
  314.    /* init global extern variables and parameters */
  315.    macheps = EPSILON; h = step; t = tol;
  316.    n = _n; x = _x; fun = _fun;
  317.    small = macheps*macheps; vsmall = small*small;
  318.    large = 1.0/small; vlarge = 1.0/vsmall;
  319.    m2 = sqrt(macheps); m4 = sqrt(m2);
  320.    ldfac = (illc ? 0.1 : 0.01);
  321.    nl = kt = 0; nf = 1; fx = (*fun)(x, n); qf1 = fx;
  322.    t2 = small + fabs(t); t = t2; dmin = small;
  323.  
  324.    if (h < 100.0*t) h = 100.0*t;
  325.    ldt = h;
  326.    for (i=0; i<n; i++) for (j=0; j<n; j++)
  327.        v[i][j] = (i == j ? 1.0 : 0.0);
  328.    d[0] = 0.0; qd0 = 0.0;
  329.    for (i=0; i<n; i++) q1[i] = x[i];
  330.    if (prin > 1) {
  331.       printf("\n------------- enter function praxis -----------\n");
  332.       printf("... current parameter settings ...\n");
  333.       printf("... scaling ... %20.10e\n", scbd);
  334.       printf("...   tol   ... %20.10e\n", t);
  335.       printf("... maxstep ... %20.10e\n", h);
  336.       printf("...   illc  ... %20u\n", illc);
  337.       printf("...   ktm   ... %20u\n", ktm);
  338.       printf("... maxfun  ... %20u\n", maxfun);
  339.    }
  340.    if (prin) print();
  341.  
  342. mloop:
  343.    sf = d[0];
  344.    s = d[0] = 0.0;
  345.  
  346.    /* minimize along first direction */
  347.    min(0, 2, &d[0], &s, fx, 0);
  348.    if (s <= 0.0)
  349.       for (i=0; i < n; i++)
  350.           v[i][0] = -v[i][0];
  351.    if ((sf <= (0.9 * d[0])) || ((0.9 * sf) >= d[0]))
  352.       for (i=1; i<n; i++)
  353.           d[i] = 0.0;
  354.    for (k=1; k<n; k++) {
  355.        for (i=0; i<n; i++)
  356.            y[i] = x[i];
  357.        sf = fx;
  358.        illc = illc || (kt > 0);
  359. next:
  360.        kl = k;
  361.        df = 0.0;
  362.        if (illc) {        /* random step to get off resolution valley */
  363.           for (i=0; i<n; i++) {
  364.               z[i] = (0.1 * ldt + t2 * pow(10.0,(double)kt)) * (random() - 0.5);
  365.               s = z[i];
  366.               for (j=0; j < n; j++)
  367.                   x[j] += s * v[j][i];
  368.         }
  369.           fx = (*fun)(x, n);
  370.           nf++;
  371.        }
  372.        /* minimize along non-conjugate directions */ 
  373.        for (k2=k; k2<n; k2++) {  
  374.            sl = fx;
  375.            s = 0.0;
  376.            min(k2, 2, &d[k2], &s, fx, 0);
  377.            if (illc) {
  378.           double szk = s + z[k2];
  379.               s = d[k2] * szk*szk;
  380.        }
  381.            else 
  382.           s = sl - fx;
  383.            if (df < s) {
  384.               df = s;
  385.               kl = k2;
  386.            }
  387.        }
  388.        if (!illc && (df < fabs(100.0 * macheps * fx))) {
  389.           illc = 1;
  390.           goto next;
  391.        }
  392.        if ((k == 1) && (prin > 1))
  393.           vecprint("\n... New Direction ...",d,n);
  394.        /* minimize along conjugate directions */ 
  395.        for (k2=0; k2<=k-1; k2++) {
  396.            s = 0.0;
  397.            min(k2, 2, &d[k2], &s, fx, 0);
  398.        }
  399.        f1 = fx;
  400.        fx = sf;
  401.        lds = 0.0;
  402.        for (i=0; i<n; i++) {
  403.            sl = x[i];
  404.            x[i] = y[i];
  405.            y[i] = sl - y[i];
  406.            sl = y[i];
  407.            lds = lds + sl*sl;
  408.        }
  409.        lds = sqrt(lds);
  410.        if (lds > small) {
  411.           for (i=kl-1; i>=k; i--) {
  412.               for (j=0; j < n; j++)
  413.                   v[j][i+1] = v[j][i];
  414.                   d[i+1] = d[i];
  415.               }
  416.               d[k] = 0.0;
  417.               for (i=0; i < n; i++)
  418.                   v[i][k] = y[i] / lds;
  419.               min(k, 4, &d[k], &lds, f1, 1);
  420.               if (lds <= 0.0) {
  421.                  lds = -lds;
  422.                  for (i=0; i<n; i++)
  423.                      v[i][k] = -v[i][k];
  424.               }
  425.        }
  426.        ldt = ldfac * ldt;
  427.        if (ldt < lds)
  428.           ldt = lds;
  429.        if (prin > 1)
  430.           print();
  431.        t2 = 0.0;
  432.        for (i=0; i<n; i++)
  433.            t2 += x[i]*x[i];
  434.        t2 = m2 * sqrt(t2) + t;
  435.        if (ldt > (0.5 * t2))
  436.           kt = 0;
  437.        else 
  438.       kt++;
  439.        if (kt > ktm)
  440.           goto fret; 
  441.    }
  442.    /*  try quadratic extrapolation in case    */
  443.    /*  we are stuck in a curved valley        */
  444.    quad();
  445.    dn = 0.0;
  446.    for (i=0; i<n; i++) {
  447.        d[i] = 1.0 / sqrt(d[i]);
  448.        if (dn < d[i])
  449.           dn = d[i];
  450.    }
  451.    if (prin > 2)
  452.       matprint("\n... New Matrix of Directions ...",v,n);
  453.    for (j=0; j<n; j++) {
  454.        s = d[j] / dn;
  455.        for (i=0; i < n; i++)
  456.            v[i][j] *= s;
  457.    }
  458.    if (scbd > 1.0) {       /* scale axis to reduce condition number */
  459.       s = vlarge;
  460.       for (i=0; i<n; i++) {
  461.           sl = 0.0;
  462.           for (j=0; j < n; j++)
  463.               sl += v[i][j]*v[i][j];
  464.           z[i] = sqrt(sl);
  465.           if (z[i] < m4)
  466.              z[i] = m4;
  467.           if (s > z[i])
  468.              s = z[i];
  469.       }
  470.       for (i=0; i<n; i++) {
  471.           sl = s / z[i];
  472.           z[i] = 1.0 / sl;
  473.           if (z[i] > scbd) {
  474.              sl = 1.0 / scbd;
  475.              z[i] = scbd;
  476.           }
  477.       }
  478.    }
  479.    for (i=1; i<n; i++)
  480.        for (j=0; j<=i-1; j++) {
  481.            s = v[i][j];
  482.            v[i][j] = v[j][i];
  483.            v[j][i] = s;
  484.        }
  485.    minfit(n, macheps, vsmall, v, d);
  486.    if (scbd > 1.0) {
  487.       for (i=0; i<n; i++) {
  488.           s = z[i];
  489.           for (j=0; j<n; j++)
  490.               v[i][j] *= s;
  491.       }
  492.       for (i=0; i<n; i++) {
  493.           s = 0.0;
  494.           for (j=0; j<n; j++)
  495.               s += v[j][i]*v[j][i];
  496.           s = sqrt(s);
  497.           d[i] *= s;
  498.           s = 1.0 / s;
  499.           for (j=0; j<n; j++)
  500.               v[j][i] *= s;
  501.       }
  502.    }
  503.    for (i=0; i<n; i++) {
  504.        if ((dn * d[i]) > large)
  505.           d[i] = vsmall;
  506.        else if ((dn * d[i]) < small)
  507.           d[i] = vlarge;
  508.        else 
  509.           d[i] = pow(dn * d[i],-2.0);
  510.    }
  511.    sort();               /* the new eigenvalues and eigenvectors */
  512.    dmin = d[n-1];
  513.    if (dmin < small)
  514.       dmin = small;
  515.    illc = (m2 * d[0]) > dmin;
  516.    if ((prin > 2) && (scbd > 1.0))
  517.       vecprint("\n... Scale Factors ...",z,n);
  518.    if (prin > 2)
  519.       vecprint("\n... Eigenvalues of A ...",d,n);
  520.    if (prin > 2)
  521.       matprint("\n... Eigenvectors of A ...",v,n);
  522.  
  523.    if ((maxfun > 0) && (nl > maxfun)) {
  524.       if (prin)
  525.      printf("\n... maximum number of function calls reached ...\n");
  526.       goto fret;
  527.    }
  528.    goto mloop;      /* back to main loop */
  529.  
  530. fret:
  531.    if (prin > 0) {
  532.          vecprint("\n... Final solution is ...", x, n);
  533.          printf("\n... ChiSq reduced to %20.10e ...\n", fx);
  534.      printf("... after %20u function calls.\n", nf);
  535.    }
  536.    
  537.    return(fx);
  538. }
  539.  
  540.